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Abstract. 

We performed kinetic simulations of diffusive shock acceleration in Type la su- 
pernova remnants (SNRs) expanding into a uniform interstellar medium (ISM). The 
preshock gas temperature is the primary parameter that governs the cosmic ray (CR) 
acceleration, while magnetic field strength and CR injection rate are secondary param- 
eters. SNRs in the hot ISM, with an injection fraction smaller than 10 are inefficient 
accelerators with less than 10 % energy getting converted to CRs. The shock structure 
is almost test-particle like and the ensuing CR spectrum can be steeper than E'^. Al- 
though the particles can be accelerated to the knee energy of lO'^^ZeV with amplified 
magnetic fields in the precursor, Alfvenic drift of scattering centers softens the source 
spectrum as steep as E^^-^ and reduces the CR acceleration efficiency. 



1. Introduction 



Most of Galactic cosmic rays up to at least the knee energy of lO^^-^eV, are believed 
to b e accelerated at SNRs within our Galaxy by diffusive shock acceleration (DSA) 
(see lHillasll2005l and references therein). In DSA theory, a small fraction of incom- 
ing thermal particles can be injected into the CR population, and accelerated to very 
high energies through their interactions with resonantly scatte ring Alfven waves in 
the converging flows across the shock ( Malkov & DrurylboOlh . Kinetic simulations 
of the CR acceleration at SNRs have shown that up to 50 % of explosion energy can 
be converted to CRs, when a fraction lO""* - 10~^ of incoming thermal particles are 



partic l es ari 

inject ed into the CR population at the gas subshock (iBerezhko & Volk:lll997l : iKangI 
I20061). This should be enough to replenish the galactic CRs escaping from our Galaxy 
with LcR ~ lO'^'erg s"'. 

Multi-band observations of nonthermal radio to y-ray emissions from several SNRs 
have been successfully explained by efficient DSA features such as high degree of shock 
compression a nd the possib le magn etic field am plification in the precursor (iReynods 



compression and the possib le magnetic neid ampimcation in the precursor (IKeynods 
Il008; Berezhk o et al.ll2009l : Morlino et all I2009h . High-resolution X-ray observations 



of several young SNRs exhibit very thin rims, indicatin g the presence of m agnetic fields 
as strong as a few 100//G downstream of the shock ( Parizot et al.ll2006i) . Moreover, 
theoretical studies have shown that efficient magnetic field amplificatio n via resonant 
and non-resonan t wave-particle interactions is an integral part of DSA dLucek & Bell 
20001 : lBelll2004l') . In the downstream region magnetic fields can be amplified by tur- 



bulence that is induced thr ough cascade of vorticity generated behind curved shocks 
(IGiacalone & Jokipiil l2007h . If there exist such amplified magnetic fields in the up- 
stream region of SNRs, CR ions with change Z might gain energies up to ^max ~ 
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lO^^-^Z eV, which may explain the all-particle CR spectrum up to the second knee at 
~ 10^^ eV with rigidity-dependent energy cutoffs. 

It has been recognized, however, that the CR spectrum at sources, A'^(£') oc 
with s <2, predicted based on the nonlinear DSA may be too flat to be consistent with 
the observed flux of CR nuclei at Earth, J{E) oc £'"2-'. Assuming an energy-dependent 
propagation path length (A oc E~^-^), a softer source sp ectrum with s ~ 2.3 - 2.4, 
is preferred by the observed data (e.g., I Ave et al1l2009h . This discrepancy could be 
reconciled if on e consider the drift of scattering centers with respect to the bulk plasma 
(ISkillinglll975l) . It reduces the velocity jump that the particles experience across the 
shock, whic h in turn softens the CR spectrum beyond the canonical test-particle slope 
(lKandl2010l : ICaprioh et al']l2010l : iPtuskin gf a/.1l20l0l) . 

Using the spherical CRASH (Cosmic-Ray Amr SHock) code, we have calculated 
the CR spectrum accelerated at SNRs from Type la SNe expanding into a uniform 
interstellar medium. 



2. Numerical Calculation 

2.1. Spherical CRASH Code in an Expanding Grid 

In the kinetic equation approach to numerical study of DSA, the following difi'usion- 
convection equation for the particle momentum distribution, f{r, p, t), is solved along 
with suitably modified gasdynamic equations 



2 t .dg 
r K{r,y)— 
or 
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where g = p^f, with f{p, r, t) the pitch angle averaged CR distribution, u^. is the wave 
speed, and y = ln(;?), and K{r,y) is the diffusion coefficient parallel to the field fines 
(ISkillin^ll97^ . 



The basic gasdynamic equations and details of the spherical CRASH code can 
be found in [Kang & Jones (2006). The advection term of Eq. ([T]) is solved by the 
wave-propagation method, as for the gasdynamic variables, except that only the entropy 
wave applies. Then the diffusion term is solved by the semi-implicit Crank-Nicholson 
scheme. In order to implement the shock tracking and AMR (Adoptive Mesh Refine- 
ment) techniques effectively in a spherical geometry, we solve the fluid and diffusion- 
convection equations in a comoving frame that expands with the outer shock. Since the 
shock is at rest and tracked accurately as a true discontinuity, we can refine the region 
around the gas subshock at an arbitrarily fine level. Moreover, the shock remains at the 
same location in the comoving grid, so the compression rate is applied consistently to 
the CR distribution at the subshock. This results in much more accurate and efficient 
low energy CR acceleration and faster numerical convergence on coarser grid spacings, 
compared to the simulations in a fixed, Eulerian grid. 

At quasi-parallel shocks, smaU anisotropy in the particle velocity distribution in 
the postshock fluid frame causes some particles in the high energy tail of the Maxw ellian 
distribution to stream upstream ( Giacalone et al. iri992l : iMaUcov & Druryll200il) . This 



thermal leakage injection process is treated numerically by adopting a phenomeno- 
logical injection scheme, in which particles above a certa in injection morn entum /jjnj 
cross the subshock and get injected to the CR population (iKang et al.ll2002h . One free 
parameter controls this process; 65 - Bo/Bj_, the ratio of the general magnetic field 
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Table 1 . Model Parameters for Supernova Remnants 

Model iiH Tq Bf, ro Ug P,, 

(cm-^^) (K) (juG) (pc) (years) (lO'^kms-') (lO^^'erg cm^^) 

WISM 0.3 3.3x10'* 30 3.19 255. 1.22 1.05 

MISM 0.03 10^ 30 6.87 549. 1.22 1.05x 10-' 

HISM 0.003 lO*- 5 14.8 1182. 1.22 1.05 xlO^^ 



WISM. eg = 0.25 : t/t„ = 0.5, 1, 3, 6, 10. 15 




Figure 1. Time evolution of the SNR model in the warm phase ISM. The CR 
distribution function at the subshock, gs = fs(p)p^, is shown as well. 



along the shock norma l, Bo, to the amplitude of postshock MHD wave turbulence, B± 
( Malkov & Volklll998h . This parameter controls the fractio n of particles inje cted into 
the CR population, ^, at the gas subshock. On the other hand, iGiacalonel (l2005h showed 
that the protons can be injected efficiently even at perpendicular shocks in fully turbu- 
lent fields due to field line meandering. 

Assuming that particles are resonantly scattered by self-generated waves, we adopt 
a Bohm-like diffusion model that represent a saturated wave spectrum, K{r, p) = Kn • 
(p/mpc){po/p), where /<-„ = nipC^ /{3eBo) = 3.13 x lO^^cm^s'^ (Bq/ jj.G)~^ . 
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WISM: £„=0.2, 0.25 MISM: £„=0.2, 0.25 HISM: e„ = 0.2, 0.25 




Figure 2. Immediate pre-subshock density, pi, post-subshock density, p2, post- 
subshock CR and gas pressure in units of the ram pressure of the unmodified Sedov- 
Taylor solution, p^U^j, and the CR injection fraction, ^ are plotted for models WISM 
(left), MISM (middle), and HISM (right). Two values of eg = 0.20 (solid Unes) and 
es - 0.25 (dashed lines), are adopted. 



2.2. Alfven Wave Transport 

The scattering by Alfven waves tends to isotropize the CR distribution in the wave 
frame, which may drift upstream at the Alfven speed, u„ = va = BqI yj4np, with re- 
spect to the bulk plasma flow, where Bq is the amplified magnetic field strength ("S killing 



19751) . In the postshock region, u„ = is assumed, since the Alfvenic turbulence in 
that region is probably relatively balanced. This Alfvenic drift reduces the velocity 
difference between upstream and downstream scattering centers compared to the bulk 
flow, leading to less efficient DSA. So the 'modified' test-particle slope can be esti- 
mated as = 3(uo - va)/(mo - va - ui), where U2 is the downstream flow speed 



(iMalkov & brurylbOOlh . Hereafter we use the subscripts '0', '1', and '2' to denote 



conditions far upstream, immediately upstream and downstream of the subshock, re- 
spectively. Thus the drift of Alfven waves in the upstream region tends to soften the 
CR spectrum from the canonical test-particle spectrum of f(p) oc p^'^, if the Alfven 
Mach number (Ma = uJva) is small. 

In addition, gas heating due to Alfven wave dissipation in the upstream region 
is considered by the term W{r, t) = -VAidPddr) where Pc is the CR pressure. This 
term is derived from a simple model in which Alfven waves are amplified by stream- 
ing CRs and dissipated locally as heat in the precursor reg ion. This effect reduces 
the subshock Mach n umber thereby reducing DSA efficiency dBerezhko & Volkll 19971 : 
Kang & Jonesll200d ). 
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Figure 3. The CR distribution at the shock, g{rs,p), the volume integrated CR 
number, Gip), and its slope, Q{p), are shown for models with eg = 0.20. 

2.3. Model Parameters for Type la Supernova Remnants 

We consider a Type la SN explosion with the ejecta mass, Mgj - 1.4Mq, expanding into 
a uniform ISM. All models have the explosion energy, Eg = 10^' ergs. The shock sonic 
Mach number is the key parameter determining the evolution and the DSA efficiency, 
while the particle injection rate and magnetic field strength (through va and k„) are 
secondary parameters. So here three phases of the ISM are considered: the warm phase 
with To = 3 X lO'^K (WISM), the intermediate phase with Tq ^ lO^K (MISM), and 
the hot phase with Tq = lO^K (HISM). The presence of amplified magnetic fields 
a few xlOOjuG downstream of shock has been indicated by very thin rims of several 
young SNRs observed in X-ray {e.g., Reynods 2008). To represent this effect we take 
the upstream field strength, Bq = 30pG for the WISM and MISM models. For the 
HISM model, however, the SNR shock is much weaker and not dominated by the CR 
pressure, so the filed amplification should be minimal. Thus the mean ISM field of 
5pG is adopted. The physical quantities are normalized by the following constants: 

ro ^ {3Mej/47:pof'^, to ^ [poti/Eof', "o ^ ro/to,Po ^ (2.34 x lO'^'^gcm"^)/!//, and 
Po = Poul- Model parameters are summarized in Table 1. Two values of eg = 0.20 and 
eg = 0.25 are adopted for inefficient and efficient injection cases, respectively. 



3. Results 

Figure 1 shows the time evolution of the WISM model with 65 = 0.25. The injection is 
very efficient with ^ w lO""*, the shock becomes dominated by Pc and the total density 
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Figure 4. Same as Figure 3 except = 0.25. 

compression peaks at pi/po « 6 at t/tg ~ 1. Afterward, the CR acceleration becomes 
less efficient as the shock slows down and becomes weaker. 

Figure 2 shows the evolution of shock properties in the different models. In the 
models with 6b = 0.2 (solid lines), the CR injection and acceleration are inefficient with 
the injected CR fraction, ^ « 10^"*, and the postshock CR pressure, PciUpqU'^t) < 0.1 
through out the entire Sedov-Taylor (ST) stage. In such inefficient acceleration regime, 
the shock remains test-particle like with P2/P0 ~ 4 and pi/po ~ 1. In the models 
with eB = 0.25 (dashed lines), on the other hand, ^ ~ 10~^, the ratio Pci/iPoU^f) 
reaches up to 0.4, and the total density compression peaks at pi/po ~ 7. These two 
ratios decrease in time, since both the sonic and Alfvenic Mach numbers decrease. The 
WISM model is the more efficient accelerator than the HISM model, because of the 
higher sonic Mach number. Despite the lower sonic Mach number, in the HISM model 
the CR acceleration is more efficient, compared to the MISM model, because of weaker 
Alfvenic drift effects due to lower Bq. 

Figures 3-4 show the CR distribution function at the shock, g{rs,p), the volume 
integrated CR spectrum, G{p) = j 47Tg{r, p)r^dr, which represents the spectrum of the 
particles confined by the shock and its slope Q{q) = -dlnG/dlnp + 4. In all mod- 
els the cutoff momentum /j^ax approaches up to ~ 10^^ - 10^^-^ eV/c. We note that, 
if the shock follows the Sedov-Taylor solution, the cutoff momentum asymptotes to 
Pmax/f^pC w 0.6\{ulto/K„) ~ 10^'^ at large t, which corresponds to Smax ~ 10^^'^ eV 
for the WISM model. Even for the inefficient injection cases with ^ < 10^^ the CR spec- 
trum at the shock exhibit concave curvatures as a consequence of momentum dependent 
diffusion across the precursor and slowing-down of the spherical shock. Softening of 
the spectrum at low energies is further enhanced by the Alfvenic drift and weakening 
of the subshock strength in time, which is more significant in the cases of stronger 
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Figure 5. Total thermal, kinetic and CR energies inside the simulation volume in 
units of the explosion energy for the models shown in Figure 2. 

Bo- However, the concavity of G{p) is much less pronounced than that of g{rs,p). In 
the inefficient injection models, the slope, Q{p), varies 4.2-4.3 at low momentum and 
3.9-4.2 at high momentum, and G{p) does not changes significantly for tjto > 6, es- 
pecially for lOp-iUpC < p < Pmax. This implies that the acceleration is nearly balanced 
by the adiabatic cooling during the late stage. So we can predict that the form of G{p) 
would remain roughly the same at much later time. In the efficient injection models 
with^> 10"^ of course, nonlinear feedback effects are more substantial. 

Figure 5 shows the integrated energies, EilE„ - An j ei{r)r^dr, where e,/,, 
and ecR are the densities of thermal, kinetic and cosmic ray energy, respectively. The 
total CR energy accelerated by t/to = 15 is Eqr/Eo = 0.15, 0.05, and 0.08 for WISM, 
MISM, and HISM models, respectively, for eg = 0.2. These are marginally consis- 
tent with the requirement that an order of 10 % of SN explosion energy needs to be 
converted to CRs to replenish the Galactic CRs. In the efficient injection models with 
es = 0.25, the CR energy fraction approaches to Eqr/Eo = 0.45, 0.25, and 0.35 for 
WISM, MISM, and HISM models, respectively. The CR acceleration in the warm ISM 
model with 65 = 0.25 may be somewhat too efficient. But one has to recall that the CR 
injection rate may depend on the mean magnetic field direction relative to the shock 
surface. In a more realistic magnetic field geometry, where a uniform ISM field is 
swept by the spherical shock, only 10-20 % of the shock surface has a quasi-parallel 
field geometry (Volk et al. 2003). Moreover, Type la SNe make up only about 15 % of 
SN explosions in the G alaxy, while core collapse SNe exploding inside wind bubbles 
may behave differendy (lHillasll2005h . 



4. Conclusion 

In general, the DSA is very efficient for strong SNR shocks, if the injection fraction 
is ^ > 10"^'^. The CR spectrum at the subshock shows a strong concavity, not only 
because the shock structure is modified nonlinearly by the dominant CR pressure, but 
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also because the spherical SNR shock slows down in time during the Sedov-Taylor 
stage. Thus the concavity of the CR spectrum in SNRs is more pronounced than that in 
plane-parallel shocks. Moreover, the Alfvenic drift in the precursor further softens the 
CR spectrum as the Alfvenic Mach number decreases. However, the volume integrated 
spectrum, G{p) (the spectrum of CRs confined by the shock including the particles in 
the upstream region) is much less concave. In the test-particle like solutions, G{p) 
approaches roughly to time-asymptotic states in the late Sedov-Taylor stage, since the 
acceleration is nearly balanced by adiabatic cooling. 

If the injection fraction is ^ > 10~^, about 25-45% of the explosion energy is 
transferred to CRs and the source CR spectrum becomes N{E) oc with s ~ 1.6-1.8 
for 10' ' < E < 10^^ e V, which might be too flat to be consistent with the observed 
CR spectrum at Earth (I Ave et al.|[2009h . If ^ < 10"'^, on the other hand, the shock 



structure is almost test-particle like with pi/po ~ 4.2 - 4.4 and the predicted source 
spectrum has a slope s « 2.0 - 2.1. However, the fraction of CR energy conversion, 
Ecr/Eq ~ 0.05 - 0.15, might be only marginally consistent with the Galactic CR 
luminosity. 

Finally, in all models considered in this study, for Bohm-like diffusion with the 
ampUfied magnetic field in the precursor, the particles could be accelerated to ^max ~ 
IQis.s^eV. The drift and dissipation of faster Alfven waves in the precursor, on the 
other hand, soften the CR spectrum and reduce the CR acceleration efficiency. 
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